library(tidyverse)
library(gWQS)
library(ggplot2)
library(cowplot)
library(knitr)
library(kableExtra)
library(gtsummary)
library(haven)
library(skimr)
library(corrplot)
library(bkmr)
library(forestplot)
library(mice)
library(rms)
library(tableone)
library(ggplot2)
library(ggcorrplot)
library(survminer)
library(sjmisc)
library(meta)
library(dplyr)
library(ggplot2)
library(survminer)
library(survival)
library(forestmodel)

## 数据提取 
# 协变量
#包括##
#年龄(连续)、性别(男性或女性)、种族或民族(非西班
#牙裔白人、非西班牙裔黑人、西班牙裔或其他)、最高教育程度成人水平(初中以下、初中、高中毕业、大学学位、
#大学毕业或以上)和收入与贫困的比率(连续)。收入/贫困比是指家庭收入与贫困指南的比值。

setwd("G:/BaiduNetdiskDownload")
# demo_c <- read_xpt("NHANES/2003-2004/Demographics/demo_c.xpt")
# demo_d <- read_xpt("NHANES/2005-2006/Demographics/demo_d.xpt")
# demo_e <- read_xpt("NHANES/2007-2008/Demographics/demo_e.xpt")
# demo_f <- read_xpt("NHANES/2009-2010/Demographics/demo_f.xpt")
# demo_g <- read_xpt("NHANES/2011-2012/Demographics/demo_g.xpt")
# demo_h <- read_xpt("NHANES/2013-2014/Demographics/demo_h.xpt")
# demo_i <- read_xpt("NHANES/2015-2016/Demographics/demo_i.xpt")
# 
# demo_all <- dplyr::bind_rows(list(demo_c,demo_d,demo_e,demo_f,demo_g,demo_h,demo_i))
# 
# #孕妇 RHD143 =1      癌症  mcq220=1排除掉  排除
# 
# 
# rhq_c <- read_xpt("NHANES/2003-2004/Questionnaire/rhq_c.xpt")
# rhq_d <- read_xpt("NHANES/2005-2006/Questionnaire/rhq_d.xpt")
# rhq_e <- read_xpt("NHANES/2007-2008/Questionnaire/rhq_e.xpt")
# rhq_f <- read_xpt("NHANES/2009-2010/Questionnaire/rhq_f.xpt")
# rhq_g <- read_xpt("NHANES/2011-2012/Questionnaire/rhq_g.xpt")
# rhq_h <- read_xpt("NHANES/2013-2014/Questionnaire/rhq_h.xpt")
# rhq_i <- read_xpt("NHANES/2015-2016/Questionnaire/rhq_i.xpt")
# 
# rhq_all <- dplyr::bind_rows(list(rhq_c,rhq_d,rhq_e,rhq_f,rhq_g,rhq_h,rhq_i))
# 
# mcq_c <- read_xpt("NHANES/2003-2004/Questionnaire/mcq_c.xpt")
# mcq_d <- read_xpt("NHANES/2005-2006/Questionnaire/mcq_d.xpt")
# mcq_e <- read_xpt("NHANES/2007-2008/Questionnaire/mcq_e.xpt")
# mcq_f <- read_xpt("NHANES/2009-2010/Questionnaire/mcq_f.xpt")
# mcq_g <- read_xpt("NHANES/2011-2012/Questionnaire/mcq_g.xpt")
# mcq_h <- read_xpt("NHANES/2013-2014/Questionnaire/mcq_h.xpt")
# mcq_i <- read_xpt("NHANES/2015-2016/Questionnaire/mcq_i.xpt")
# 
# mcq_all <- dplyr::bind_rows(list(mcq_c,mcq_d,mcq_e,mcq_f,mcq_g,mcq_h,mcq_i))
# 
# 
# # ，在血清中测量了四种PFAS化合物包括
# #全氟辛酸(PFOA)PFAS lbxnfoa、  只在 PFAS有.   在PFC l24pfc 叫 LBXPFOA  perfluorooctanoic acid
# #全氟辛烷磺酸(PFOS  lbxnfos)、 只在 PFAS有. 在PFC l24pfc 叫 LBXPFOS  perfluorooctane sulfonic acid
# #全氟己烷磺酸(PFHxS lbxpfhs)和 只在 PFAS 和 PFC 和 l24pfc 有  perfluorohexane sulfonic acid
# #全氟壬酸(PFNA lbxpfna)。 全有 PFAS 和 PFC 有  perfluorononanoic acid
# # 两部分获取 部分数据存在不同模块里面
# # 
# # L24PFC  2003-2004
# #PFC 2005-2012
# #PFas 2013-2016
# # 注：使用2013-2014年NHANES的剩余血清标本测量了另外四种全氟烷基和多氟烷基物质（PFAS）
# #，即PFOA和PFOS的四种异构体。这些异构体数据应与本数据集（PFAS_H）结合使用，并可在本网站的全氟烷基和多氟烷基物质（PFAS）-直链和支链全氟辛烷磺酸和PFOA异构体（SSPFAS_H）数据集中找到。
# #从NHANES 2015-2016开始，这些异构体将包含在全氟烷基和多氟烷基物质（PFAS）数据集中。
# #您也可以使用以下链接访问这些异构体的数据集： 
# # 2013-2014年PFAS（以前称为PFC）数据集中计算出的异构体总和与NHANES前几个周期报告的总水平相当。
# pfas_c <- read_xpt("NHANES/2003-2004/Laboratory/l24pfc_c.xpt")
# pfas_d <- read_xpt("NHANES/2005-2006/Laboratory/pfc_d.xpt")
# pfas_e <- read_xpt("NHANES/2007-2008/Laboratory/pfc_e.xpt")
# pfas_f <- read_xpt("NHANES/2009-2010/Laboratory/pfc_f.xpt")
# pfas_g <- read_xpt("NHANES/2011-2012/Laboratory/pfc_g.xpt")
# # 这里面数据是
# # linear PFOA n-PFOA -> LBXNFOA
# # sum of branched isomers of PFOA  Sb-PFOA -> LBXBFOA
# # linear PFOS n-PFOS ->  LBXNFOS
# # the sum of monomethyl branched isomers of PFOS ->
# pfas_h <- read_xpt("NHANES/2013-2014/Laboratory/pfas_h.xpt")
# pfas_i <- read_xpt("NHANES/2015-2016/Laboratory/pfas_i.xpt")
# 
# # 进行变量添加转换后面用到 2023年7月24日10:41:46 搞不定,13年和15年数据不一致 变量也不对应 无法搞 这样
# 
# pfas_all <- dplyr::bind_rows(list(pfas_c,pfas_d,pfas_e,pfas_f,pfas_g,pfas_h,pfas_i))
# 
# # BMI 
# 
# bmx_c <- read_xpt("NHANES/2003-2004/Examination/bmx_c.xpt")
# bmi_d <- read_xpt("NHANES/2005-2006/Examination/bmx_d.xpt")
# bmi_e <- read_xpt("NHANES/2007-2008/Examination/bmx_e.xpt")
# bmi_f <- read_xpt("NHANES/2009-2010/Examination/bmx_f.xpt")
# bmi_g <- read_xpt("NHANES/2011-2012/Examination/bmx_g.xpt")
# bmi_h <- read_xpt("NHANES/2013-2014/Examination/bmx_h.xpt")
# bmi_i <- read_xpt("NHANES/2015-2016/Examination/bmx_i.xpt")
# bmi_all <- dplyr::bind_rows(list(bmx_c,bmi_d,bmi_e,bmi_f,bmi_g,bmi_h,bmi_i))
# 
# # 饮食中的叶酸摄入
# # 第1天
# dr1tot_c <- read_xpt('NHANES/2003-2004/Dietary/dr1tot_c.xpt')
# dr1tot_d <- read_xpt('NHANES/2005-2006/Dietary/dr1tot_d.xpt')
# dr1tot_e <- read_xpt('NHANES/2007-2008/Dietary/dr1tot_e.xpt')
# dr1tot_f <- read_xpt('NHANES/2009-2010/Dietary/dr1tot_f.xpt')
# dr1tot_g <- read_xpt('NHANES/2011-2012/Dietary/dr1tot_g.xpt')
# dr1tot_h <- read_xpt('NHANES/2013-2014/Dietary/dr1tot_h.xpt')
# dr1tot_i <- read_xpt('NHANES/2015-2016/Dietary/dr1tot_i.xpt')
# 
# dr1tot.data.file <- dplyr::bind_rows(list(dr1tot_c, dr1tot_d,dr1tot_e,dr1tot_f,dr1tot_g,dr1tot_h,dr1tot_i))
# dr1tot.data <- dr1tot.data.file[,c('SEQN', 'DR1TFA')]
# 
# # 第2天
# dr2tot_c <- read_xpt('NHANES/2003-2004/Dietary/dr2tot_c.xpt')
# dr2tot_d <- read_xpt('NHANES/2005-2006/Dietary/dr2tot_d.xpt')
# dr2tot_e <- read_xpt('NHANES/2007-2008/Dietary/dr2tot_e.xpt')
# dr2tot_f <- read_xpt('NHANES/2009-2010/Dietary/dr2tot_f.xpt')
# dr2tot_g <- read_xpt('NHANES/2011-2012/Dietary/dr2tot_g.xpt')
# dr2tot_h <- read_xpt('NHANES/2013-2014/Dietary/dr2tot_h.xpt')
# dr2tot_i <- read_xpt('NHANES/2015-2016/Dietary/dr2tot_i.xpt')
# 
# dr2tot.data.file <- dplyr::bind_rows(list(dr2tot_c, dr2tot_d,dr2tot_e,dr2tot_f,dr2tot_g,dr2tot_h,dr2tot_i))
# dr2tot.data <- dr2tot.data.file[,c('SEQN', 'DR2TFA')]
# 
# # 合并第1天&第2天的变量
# dr.data <- merge(dr1tot.data, dr2tot.data)
# # 进行计算平均值 
# # day1、day2 都有，则取2天平均，否则，取第一天的值
# # 先计算2天的平均值，其中第二天为 NA 的，平均值也是 NA
# totaldrFA <- apply(dr.data[,c('DR1TFA', 'DR2TFA')], 1, mean)
# dr.data$totaldrFA <- totaldrFA
# 
# # 把第二天为 NA 的值计算出来，用第一天的值作为平均值
# day.2.na.index <- which(is.na(dr.data$totaldrFA))
# totaldrFA[day.2.na.index] <- dr.data$DR1TFA[day.2.na.index]
# # 添加到paper.data中，再次确认
# dr.data$totaldrFA <- totaldrFA
# 
# # 结束计算平均值
# 
# 
# # 通过文章里面描述 是多个不同周期需要合并或者字段不统一 需要进行合并 麻烦了 ng/ml
# # 血清叶酸
# # FOLFMS	lbdfotsi
# #2003 2006 L06NB_C d 数据 2023年7月25日16:06:52 
# # 血清叶酸 serumfo  红细胞叶酸 rbcfo 这是我自定义的 部分不同字段都合并到这个字段来 统一处理
# 
# folfms_c <- read_xpt("NHANES/2003-2004/Laboratory/l06nb_c.xpt")
# # 统一变量
# folfms_c$serumfo <- folfms_c$LBXFOL
# folfms_c$rbcfo <- folfms_c$LBXRBF
# folfms_d <- read_xpt("NHANES/2005-2006/Laboratory/folate_d.xpt")
# folfms_d$serumfo <- folfms_d$LBXFOL
# folfms_d$rbcfo <- folfms_d$LBXRBF
# 
# folfms_e <- read_xpt("NHANES/2007-2008/Laboratory/folate_e.xpt")
# folfms_e$serumfo <- folfms_e$LBDFOL
# folfms_e$rbcfo <- folfms_e$LBDRBF
# 
# folfms_f <- read_xpt("NHANES/2009-2010/Laboratory/folate_f.xpt")
# folfms_f$serumfo <- folfms_f$LBDFOL
# folfms_f$rbcfo <- folfms_f$LBDRBF
# 
# folfms_g <- read_xpt("NHANES/2011-2012/Laboratory/folfms_g.xpt")
# folfms_g$serumfo <- folfms_g$LBDFOT
# folfms_h <- read_xpt("NHANES/2013-2014/Laboratory/folfms_h.xpt")
# folfms_h$serumfo <- folfms_h$LBDFOT
# folfms_i <- read_xpt("NHANES/2015-2016/Laboratory/folfms_i.xpt")
# folfms_i$serumfo <- folfms_i$LBDFOT
# 
# folfms_all <- dplyr::bind_rows(list(folfms_c,folfms_d,folfms_e,folfms_f,folfms_g,folfms_h,folfms_i))
# # 红细胞叶酸 上面部分拿到数据了这里就只取2011-2016的红细胞叶酸
# #FOLATE  LBDRFO 
# folate_g <- read_xpt("NHANES/2011-2012/Laboratory/folate_g.xpt")
# folate_g$rbcfo <- folate_g$LBDRFO
# folate_h <- read_xpt("NHANES/2013-2014/Laboratory/folate_h.xpt")
# folate_h$rbcfo <- folate_h$LBDRFO
# folate_i <- read_xpt("NHANES/2015-2016/Laboratory/folate_i.xpt")
# folate_i$rbcfo <- folate_i$LBDRFO
# folate_all <- dplyr::bind_rows(list(folate_g,folate_h,folate_i))
# 
# output <- plyr::join_all(list(demo_all,rhq_all,mcq_all, pfas_all,folfms_all,bmi_all,folate_all,dr.data),
#                          by='SEQN',type='full')
# 
# 
# 
# 
# # 提取数据 删除没用数据
# # 我们纳入了年龄在12岁及以上的调查参与者他们拥有血清中PFAS浓度、
# #红细胞和血清中叶酸生物标志物以及协变量(如协变量部分所述)的完整数据
# # 12岁以上的
# output12 <- output[which(output$RIDAGEYR >= 12), ]
# 
# 
# 
# # output12 <- output
# # 这块数据不对 始终对不起来
# #PFAS 完整数据的  四种数据衍生物进行删除na的
# #全氟辛酸(PFOA)PFAS lbxnfoa、 只在 PFAS有 在PFC  叫 LBXPFOA  perfluorooctanoic acid
# #全氟辛烷磺酸(PFOS  lbxnfos)、 只在 PFAS有 在PFC 叫 LBXPFOS  perfluorooctane sulfonic acid
# #全氟己烷磺酸(PFHxS LBXPFHS)和 只在 PFAS 和 PFC有  perfluorohexane sulfonic acid
# #全氟壬酸(PFNA LBXPFNA)。 全有 PFAS 和 PFC L24PFC有  perfluorononanoic acid
# output12 <- output12[which(!is.na(output12$LBXPFOA)),]
# output12 <- output12[which(!is.na(output12$LBXPFOS)),]
# output12 <- output12[which(!is.na(output12$LBXPFHS)),]
# output12 <- output12[which(!is.na(output12$LBXPFNA)),]
# #暂且这样吧 其他不知道为什么对不起来
# # table(output12$RIDAGEYR)
# 
# 
# 
# # 删除不含有血清叶酸和红细胞叶酸  的数据 
# # 血清叶酸是五个维度累加得到的 采用我的自定义变量  血清叶酸 serumfo  红细胞叶酸 rbcfo
# 
# output12 <- output12[which(!is.na(output12$serumfo)),]
# output12 <- output12[which(!is.na(output12$rbcfo)),]
# 
# # 删除缺失协变量的 年龄, sex, race, BMI, PIR, 教育程度 
# output12 <- output12[which(!is.na(output12$RIDAGEYR)),]
# output12 <- output12[which(!is.na(output12$RIAGENDR)),]
# output12 <- output12[which(!is.na(output12$RIDRETH1)),]
# output12 <- output12[which(!is.na(output12$BMXBMI)),]
# output12 <- output12[which(!is.na(output12$INDFMPIR)),]
# 
# output12lower <- output12[which(output12$RIDAGEYR < 20),]
# 
# output20 <- output12[which(output12$RIDAGEYR >= 20), ]
# 
# output1220edu <- output20[which(!is.na(output20$DMDEDUC2)),]
# 
# output12<-dplyr::bind_rows(list(output12lower,output1220edu))
# 
# # output12back <- output12
# 
# # 排除孕妇 曾经得过癌症的排除
# 
# #孕妇 RHD143 =1 癌症  mcq220=1
# 
# output12 <- output12[-which(output12$RHD143==1),]
# output12 <- output12[-which(output12$MCQ220==1),]